


1. Report No. 2. Government Accession No. 

NASA TR R-436 

3. Recipient's Catalog No. 

4. Title and Subtitle 

A CUBIC SPLINE APPROXIMATION FOR PROBLEMS 
IN FLUID MECHANICS 

5. Report Date 

October 1975 

6. Performing Organization Code 

7. Author(s) 

Stanley G. Rubin and Randolph A. Graves, Jr. 

8. Performing Organization Report No. 

L-9929 

10. Work Unit No. 

506-26-20-05 

9. Performing Organization Name and Address 

NASA Langley Research Center 
Hampton, Va, 23665 

11. Contract or Grant No. 

13. Type of Report and Period Covered 

Technical Report 

12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Washington, D.C. 20546 

14. Sponsoring Agency Code 


15. Supplementary Notes 


Stanley G. Rubin is professor of aerospace engineering, Polytechnic Institute of New York, 
Farmingdale, N.Y. Work performed as visiting professor at Old Dominion University, Norfolk, 
Va. 


16. Abstract 

A cubic Spline approximation is presented which is suited for many fluid-mechanics 
problems. This procedure provides a high degree of accuracy, even with a nonuniform 
mesh, and leads to an accurate treatment of derivative boundary conditions. The trun- 
cation errors and stability limitations of several implicit and explicit integration schemes 
are presented. For two-dimensional flows, a spline -alternating -direction-implicit (SADI) 
method is evaluated. The spline procedure is assessed, and results are presented for the 
one -dimensional nonlinear Burgers^ equation, as well as the two-dimensional diffusion 
equation and the vorticity-stream function system describing the viscous flow in a driven 
cavity. Comparisons are made with analytic solutions for the first two problems and with 
finite -difference calculations for the cavity flow. 


17. Key Words (Suggested by Author(s)| 

Cubic splines 
Burgers’ equation 
Cavity flow 
Numerical analysis 
Splines under tension 


18. Distribution Statement 

Unclassified - Unlimited 


Subject Category 12 


19. Security Oassif. (of this report) 
Unclassified 


20. Security Classif. (of this page) 

21. No. of Pages 

Unclassified 

93 


22, Price* 

$4.75 


For sale by the National Technical Information Service, Springfield, Virginia 221 61 





CONTENTS 


SUMMARY 

INTRODUCTION 

SYMBOLS 

SPLINE FORMULATION 

Basic Spline Theory 

Splines for Solving Partial Differential Equations .... 

Linear Burgers’ equation 

Nonlinear Burgers' equation 

Two-dimensional equation 

TRUNCATION ERROR 

Theory for Cubic ^lines 

Examples of Truncation Error Using Burgers' Equation 

Second -order spatial derivative 

First-order temporal derivative 

Diffusion only 

Convection only 

Complete equation 

STABILITY 

General Development for Linearized Burgers' Equation . 

Explicit convection and diffusion 

Explicit convection and implicit diffusion 

Implicit convection and diffusion 

Development for SADI Procedure 

SPLINE CURVE FITTING 

Resolution 

Splines Under Tension 

Divergence Form 

BURGERS' EQUATION 

Numerical Procedure 

Implicit method 

Two-step method 

Discussion of Numerical Results . 

DIFFUSION EQUATION 


Page 

1 

1 

3 

6 

6 

8 

8 

10 

11 

12 

12 

13 

13 

15 

15 

15 

16 

17 

17 

18 
19 
21 
21 

23 

23 

26 

28 

29 

29 

29 

30 
30 

32 


iii 


Page 

INCOMPRESSIBLE FLOW IN A CAVITY 33 

Numerical Procedure 33 

Stream function 33 

Vorticity 33 

Discussipn of Numerical Results 35 

CONCLUDING REMARKS 37 i 

APPENDIX - KREISS FOURTH-ORDER FINITE-DIFFERENCE METHOD 39 

REFERENCES 40 

TABLES 42 

FIGURES 79 


iv 


I li liill lllii lii. illL: 



A CUBIC SPLINE APPROXIMATION FOR PROBLEMS 
IN FLUID MECHANICS 

Stanley G. Rubin* and Randolph A. Graves, Jr. 

Langley Research Center 

SUMMARY 

A cubic spline approximation is presented which is suited for many fluid-mechanics 
problems. This procedure provides a high degree of accuracy, even with a nonuniform 
mesh, and leads to an accurate treatment of derivative boundary conditions. The trun- 
cation errors and stability limitations of several implicit and explicit integration schemes 
are presented. For two-dimensional flows, a spline -alternating-direction-implicit (SADI) 
method is evaluated. The spline procedure is assessed, and results are presented for the 
one -dimensional nonlinear Burgers' equation, as well as the two-dimensional diffusion 
equation and the vorticity-stream function system describing the viscous flow in a driven 
cavity. Comparisons are made with analytic solutions for the first two problems and 
with finite -difference calculations for the cavity flow. 

INTRODUCTION 

The numerical treatment of many problems in fluid mechanics is complicated by 
three conditions: (1) local singular regions where the flow gradients are much larger 

than typically found over the remainder of the domain, e.g., in the limit of large Reynolds 
number (particular examples of this singular behavior are given by shock waves, bound- 
ary and shear layers, entropy layers, etc.); (2) curvilinear boundaries that do not pass 
directly through the nodal points of a fixed uniform mesh. (Here, both geometric sur- 
faces, as well as discrete shock waves, are referred to as they appear in a numerical 
shock fitting procedure with a fixed mesh and moving shock, Moretti, ref. 1), and 
(3) derivative boundary conditions as occur for vorticity or pressure (see Roache, ref, 2). 

In some cases, coordinate transformations can alleviate the difficulties associated 
with conditions (1) and (2), but this generally requires a priori knowledge of the local or 
asymptotic flow behavior, which is not always available. Moreover, suitable transfor- 
mations are difficult to formulate if multiple shocks or other singular regions appear, if 

♦Professor of aerospace engineering. Polytechnic Institute of New York, 
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a number of geometric configurations must be considered simultaneously, or if a singular 
region is multilayered, i.e., singular regions within singular regions. Some examples of 
the latter are the trailing-edge boundary layer (Messiter, ref. 3), the laminar sublayer 
within a turbulent boundary layer (ref. 4), oscillating boundary layers (Ackerberg, ref. 5), 
and corner boundary regions (Rubin, ref. 6). 

The accuracy of a numerical calculation can be improved by suitable mesh reduc- 
tion or by increasing the order of the truncation error. However, higher order methods 
generally require the introduction of additional nodal points in the discretization formu- 
las, thereby increasing the coupling in the system of algebraic difference equations. For 
implicit methods the number of nonzero entries in the Inversion matrix is increased so 
that the tridiagonal form associated with a three-point formulation no longer occurs. 

Since the very efficient tridiagonal inversion algorithms can no longer be applied, a sig- 
nificant increase in computer time results. 

Higher order discretizations can also be used for accurately representing deriva- 
tive boundary conditions (Briley, ref. 7). However, these may be inadequate if the mesh 
dimension is too large; i.e,, if the local surface gradients are o(a- 1) and the mesh 
dimension is 0(A), the accuracy will remain poor regardless of the number of terms 
retained in a Taylor series expansion. In many problems a surface layer grows from 
zero thickness Initially to some finite thickness in a steady state. In the initial stages, 
inaccuracy near the boundary can lead to a divergence that is suppressed only by an 
under -relaxation procedure (Bozeman and Dalton, ref. 8). 

Uniform mesh reduction improves accuracy but results in a significant increase in 
the number of algebraic difference equations and is particularly inefficient from the point 
of view of computer storage and calculation time. A nonuniform mesh that is adjusted to 
reflect the appearance of singular regions and irregular boundaries should be optimal. 
Unfortunately, with a three -point finite -difference approximation the order of the trun- 
cation error will be significantly decreased with even a moderate variation in the mesh 
dimension (Crowder and Dalton, ref. 9), Therefore, the expected increase in accuracy 
associated with mesh reduction is not achieved. 

The present paper describes a cubic spline procedure for the solution of second- 
order quasi -linear partial differential equations in one or two spatial dimensions. A 
finite -difference discretization is used for the marching or time-like direction. Unlike 
a finite-element or Galerkin procedure, there are no quadratures to evaluate, and the 
coefficient matrix is tridiagonal. Implicit and explicit spline fitting is examined for a 
one -dimensional Burgers’ equation; a spline-alternating-direction-implicit (SADI) pro- 
cedure is formulated for two-dimensional flows. The spline approximation is second- 
order accurate, even with relatively large variations in the mesh, so that singular regions 
and irregular boundaries can be considered without loss of accuracy and with a minimum 
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of computer storage and time. Moreover, for invlscid flows where the system of differ- 
ential equations becomes first order, the spline procedure is third-order accurate with 
a nonuniform mesh and of fourth-order accuracy with a uniform mesh. This result is 
consistent with the increased accuracy of lower order derivatives in a spline curve fit 
(Ahlberg, Nilson, and Walsh in ref. 10). For a uniform mesh, a particular combination 
of spline’s and finite differences results in a fourth-order accurate procedure for viscous 
flows as well. The tridiagonal form is maintained. 

Since the spline approximation provides a direct relation between the derivatives 
and the functional values evaluated at the nodal points, a finite -difference discretization 
is unnecessary. Derivative boundary conditions are imposed directly without incurring 
large local discretization errors due to inaccurate higher order one-sided difference 
approximations. This represents a significant advantage of the spline technique over 
conventional finite-difference procedures. Finally, unlike finite -difference or Galerkin 
techniques, with a spline approximation there appears to be no particular advantage 
gained with the divergence form of the equations. This fact, previously noted by Douglas 
and Dupont in reference 11 in their collocation procedure, can prove extremely important 
for flow problems where shock waves are captured during the numerical computation. 

The spline formulation and the procedure for solving second-order quasi -linear 
partial differential equations are reviewed; the truncation errors for second derivatives, 
i.e., diffusion, and first derivatives, i.e., convection, are explicitly elucidated; the sta- 
bility of explicit, implicit, and combined two-step procedures, with a uniform mesh, is 
discussed for a linearized Burgers’ equation in one dimension and with the SADI proce- 
dure in two dimensions; and the concepts of splines under tension as a smoothing proce- 
dure are reviewed. The effects of tension, mesh variation, and divergence form on the 
resolution of a spline curve fit are discussed. Also, results are presented for the one- 
dimensional nonlinear Burgers’ equation, the two-dimensional diffusion equation, and the 
viscous flow in a driven cavity. Comparisons are made with exact solutions and finite- 
difference solutions where available. 


SYMBOLS 

end points of interval or dimensions of cavity 
ai,bi,ci,di scalar coefficients in equation (10) 
c stability coefficient, u At/h 

65 truncation error 
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g(x) 

hij,hi,h 

h 

fi,mi 

Li,Mi,P 

frii 

R 

Rc 

t 

At 

Ti,Pi 

Tr,Pr 

u,v 

u 

U 

Vi 

x,y,z 

Xi 

Y,Z 


initial conditions on velocity 
mesh width, xj - Xi_i 
average mesh width 
mesh width, yj - yi_i 

spline first derivative in y- and x-direction, respectively 

spline second derivative in y-, x-, and z -direction, respectively 

spline first derivative of u2/2 

Reynolds number 

cell Reynolds number, \ih/p 

time 

time step increment 

matrices in equation (27) 

matrices in SADI stability analysis 

velocity in x- and y-direction, respectively 

coefficient in linear Burgers' equation 

wave speed in Burgers' equation 

vector 

spatial coordinates 
nodal points 
normalized coordinates 
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coefficient in spline solution procedure 


Ai 

^ vorticity 

77 transformed coordinate, x - ut 

0,01,02 finite -difference scheme (0 for explicit, 1/2 for Crank -Nicolson, 1 for 
implicit) 

eigenvalue 

1/ kinematic viscosity 

(T tension factor 

crj spacing factor 

T fictitious time 

At fictitious time step increment 

u amplification factor (see eq. ( 27 )) 

\j/ stream function 

u) wave number 

Superscripts; 

n time step number 

s fictitious time step number 

Subscripts: 

indices for x- and y-direction, respectively 
number of nodes on [a,l^ excluding the boundaries 
differentiation with respect to time 
differentiation with respect to spatial coordinates 


N 

t 


SPLINE FORMULATION 


Basic Spline Theory 

Consider a mesh with nodal points (knots) xj such that 
a = xo < XI < X 2 • • ■ < x^ < xjg+i = b 

and with 

hj = xj - Xi_i > 0 

Consider a function u(x) such that at the mesh points xj 
u(xi) = Ui 


The cubic spline is a function S^(u;x) = S^(x) which is continuous together with its 
first and second derivatives on the interval [a,^, corresponds to a cubic polynomial in 
each subinterval Xj_i S x 5 X[, and satisfies S^(ui;xi) = up 

If the function u(x) belongs to C^ja,!^, it has been shown that the spline function 
Sa(x) approximates u(x) at all points in |a,l^ to fourth order in maximum hp First 
and second derivatives of S^(x) approximate u'(x) and u"(x) to third and second 
order, respectively. See Ahlberg, Nilson, and Walsh in reference 10 for detailed proofs 
of convergence and for a discussion concerning the relationship of this spline approxi- 
mation with a mechanical spline. 


If S^(x) is cubic on [xi_i,x^, then 


/X /Xi - x\ /X - X: 1 

S’AW = + Mil— ^ 


V hi 


where Mj = S'^(xj). 


Integrating twice leads to the useful interpolation formula on [xi_i,x^ as follows: 

C / \ Tv/r (^i - (x - Xi_i)3 ( Mi_ihf\xj-x ( MihfV-xs i 

* M, ^ * (u, - (1) 


The constants of integration have been evaluated from S^(xi) = u[ and S^(xi_j) = u[_i 
where S^(x) on (i,i+^ is obtained with i+1 replacing i in equation (1). 

The unknown derivatives Mj are related by enforcing the continuity condition on 
S^(x). With S’(xj') = mf on (1-1, i] and S)^(xf) = mj" on [i,i+^, 

mf = mf = m[ 
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For i = 1, . . N, 


Mi 1 + ^ 

6 


+ ^i+1 


Mi 




^i - ^i-l 

hi 


( 2 ) 


Additional relationships obtained from equations (1) and (2), which prove useful later 
herein, are listed as follows: 


i-i. 

hi 


1 + 


2fi 


hi+1, 


mi + 


hi+l 


mi+i = 


3(ui+l-ui) 3(ui-ui_i) 


h?+l 


ht 


(3) 


^i+i - mi 


hj+1 

2 


(Mi + Mi+i) 


(4) 


hi _ _ hi U[ - u 

"■i = T + y “>-i ~ 


Ui-1 


or 


mi = - 


hi+l 


Mi 


hi+l 


Mi+i + 


Uj+1 - ^i 
hi+l 


(5) 


( 6 ) 


Therefore, given the values ui, the equations (2) and (3) with appropriate boundary 
conditions form a closed system for mi and Mi; and with equation (1) the values 
Sa(x) can be found at all intermediate locations. Equation (2) or (3) lead to a system 
of N equations for the N + 2 unknowns Mi or mi, respectively. The additional 
two equations are obtained from boundary conditions on mQ and m^+l or Mq and 
Mn+ 1- The resulting tridiagonal system for Mi or mi is diagonally dominant and 
solved by an efficient inversion algorithm (see Ahlberg, Nilson, and Walsh, ref. 10 or 
Keller, ref. 12). Note that if Mq and M^^.! are given so that all Mi (i = 1, . , ., N) 
are determined from equation (2), then mg and m^j^j are found from equation (5) 
or (6). If mo and m^+i or Amg + BMq and Cm^+i + are prescribed, 

then and mg are eliminated with equations (5) and (6) in favor of Mj,j and 

Mjvf+i and Mq and Mj, respectively. This gives a relation of the form 

EMq + FMj = G(uo,ui) 


and 

HMn + = K(un,un+i) 

where A to F, H, and J are constants and G and K are functions of the velocity 
u. These two conditions with equation (2) then close the system. 
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Splines for Solving Partial Differential Equations 

If the values ui are not prescribed but represent the solution of a quasi -linear 
second-order partial differential equation 

ut = f(u,Ux,Uxx) 

then an approximate solution for u^ can be obtained by considering the solution of 

(ut)j = f(ui,mi,Mi) 

where the time derivative is discretized in the usual finite -difference fashion: 


uP+1 - uP 

1 


At 


i= (1 - 0)f" + 0fn+l 


Linear Burgers' equation . - Consider the linear Burgers' equation 
ut + UUx = t/Uxx 

where u = u(x,t) and v = !^(x,t). Therefore, 

uP+1 = uP - At (l - 0i)uPmP + ejuP+lmP+lJ + At^(l - 02)*^pMP + 

With equations (2) and (3) a system of 3N equations for 3(N + 2) unknowns is 
obtained; the system can be written as 


AiVP+jl + B^VP+I + CiVP+/ = DjVP 


where 


Ai = 


0 0 0 

0 ^ 
hi 6 

O 1 

0 


hf 


h: 


(7) 


(8) 


(9a) 


(9b) 


Bi = 


«0 


1 + 1 


^i ^i+1 


T+1 


«1 

0 


^+1 ^i 


«2 

hj + hj+i 


3 

0 


(9c) 
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0 


0 


hi+1 


-3 

1 

_h?+l 

hi+1 

Pq P\ 

P2 

0 0 

0 

_0 0 

0 


Vi = [ui,mi,M^'r 

QfO = 1 

«1 = At 

®2 ~ 

P0= 1 

pi = -(1 -®iK 

P2 = (l - At 


(9d) 


(9e) 


(9f) 


(9g) 


Initial conditions are specified such that u(x,0) - g(x). If boundary conditions are 
specified on u(a,t) = ri(t) and u(b,t) = rjW, then ug+' and uRVi are given as rj(t) 
and r 2 (t), respectively. With derivative boundary conditions ux(a,t) = sj(t) or 
ux(b,t) = S 2 (t), mg and are prescribed as s^{t) and S 2 (t), respectively. 

From equation (8), ug+1 is given as a function of M"+l and mg+l and u"+ is 
given as a function of m^+l and M"+1; from equation (4), m^+l is given as a func- 
tion of Mg+1, M^+^and mg+1. 

With these relations for u^+^ and m^+1, and with either ug+^ or mg+^ spec- 
ified, equation (5) or (6) provides a linear relationship between Mg and M^. A similar 

result can be obtained for and The system is now closed and system (9) 

can be solved by the tridiagonal algorithm previously noted. An analogous procedure 
determines the appropriate relationships between ug and uj and ujg and uj^^l for 
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mg and specified, or mg and mj and mj^ and for ug and 

specified. 


The system (9) can be reduced by substitution of Ui and mj as functions of 
into a single tridiagonal system for Mi. The resulting equations for Mi 
(i = 1, , , ,, N) are, with 9 = 9i = 82 , 


„ TiTn+1 

aiMi.i 


+ biMf+^ 


+ CiMf;/ = di 


where 


(10) 


ai 


= h . 9(hilhj:i + 

6 ' y 3Ai hiAiy 


c. = ^ + 9(^h±ll^ . _n±i_X^^ 

^6 y 3Ai+i hi+iAi+jy 


_ hj + hj^.i n/^i+1 + 26i 26i + 6i_i 




3Ai+i 


vn+1 


3Ai ^i+lbi+i '^ihi/ 


,, /u_ i+i - Uj Uj^ Ui-l\*^ /, px/26imi - 26i+imi+i 26imi - 26i_imi_i 

‘ W+l^i+1 hiAi ; hi+iAi+i + h[Al 


n+lMj+i - yjMj yjMj - yj_iMi_i\a 


Vl^i+l 


hiAi 


where 26i = ui At, ri = At, and Ai = 1 + 2hf^(6i - 6i_i)0. 

The boundary conditions for Mg and Mj^^j are obtained in the same manner as 
outlined previously. A tridiagonal relationship similar to equation (10) can also be found 
for m""^ or uf"*"^, although the manipulation is somewhat more tedious. 

Nonlinear Burgers* equation . - If the governing equation is nonlinear 

ut + UUx = I^Uxx 

then the spline formulations gives, with 0 = 0j = 02 , 

= up + 0 At(-u^m^ + + (1 - 0) At^-u^mj + yjMj)” (11) 

If quasi-linearization is used for (uimi)"'*'^, it is found that 
(uj[mi)>^+l = uPmP'*'^ + uP'*'^mp - upmp 
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and therefore 

uf+Hi * 9 - ■'r’wr^) + - 9)9" iw" - » - 


( 12 ) 


With equation (12) in place of equation (8), the system (9) is of the same form but with the 
following modifications: 


o-Q = 1 + 0 AtmP 


(13) 


ai = 0uP At ^ 

pj = -(1 - 20)uP At 

Two-dimensional equation. - For equations with two space dimensions such that 
Ut = f(u,Ux,Uy,Uxx>Uyy) 

a spline-alternating-direction-implicit (SADI) formulation is developed. The two-step 
procedure, with quasi -linearization or some other iterative process used for nonlinear 
terms, is of the following form: For step 1, 




"1] 

= Uij +y f(uij 2,^ij,L?j^ 

) 


and for step 2, 






u9.+l 

»3 

_/4^At 



where 

^ii 

and 

Lj 4 are the spline approximations to ^Uy^ . 

.. and (uyy).., 1 

Therefore, 

in two dimensions with hjj = 

and kj 

ij = Yij - yi,j-D 

hfjVi 

-i,i ^ 

2(hfjl 

‘ + = 3hr2, j 

("i+i,) 

- U|jj + 3h|j2^Uij 

hijMi. 

i,j + 

2(hij . 

^ + hi+ijMi+ij 1= (ui+ij 

Uij) - 6hr.l(uij - 

and 






k.-.lf. 
i] 1,: 

i-i"" 

2(k,f 


i,j+l ■ 

ujj) + 3ky2^ujj - 


(14a) 


(14b) 


(15a) 


(15b) 


(16a) 
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+ 2(kjj + - 6kj-j_^i(Uij^l - u^) - Sky^^uy 


i,J 


■l) 


(16b) 

with expressions similar to equation (4) to equation (6) relating to Mjj and H-- 

to Lij. '' 

If cross derivatives such as Uxy appear in the governing system, the spline 
approximation for these terms is found from equation (16a), with mij replacing uy 
and fjj replacing The solutions fjj are the necessary spline fits to Uxy. 
Alternatively one could replace u^j with Cy and my with mij in equation (15a). 
The result ?ij = mij should be correct to third order in maximum (hjj, kij). 


TRUNCATION ERROR 


Theory for Cubic Splines 

For interior points, the spatial accuracy of the spline approximation can be 
directly estimated from the formulas (2) and (3) or the equivalent two-dimensional rela- 
tionships (15a) and (15b). Expanding mjj, Mij, and uy in Taylor series and assum- 
ing the necessary continuity of derivatives for u(x,y) gives 



Fyfe (ref. 13) has presented similar relations, for constant h^, in his collocation analy- 
sis of cubic splines for the solution of two-point boundary value problems. 

Therefore, the spline approximation with nonuniform mesh is second-order accu- 
rate for Mjj and third-order for m^j. For a uniform mesh, mjj is fourth-order 
accurate; and with h^j = h, 

2 

Mij = (uxx)ij - (uxxxx)ij y 2 + (18) 
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The standard three -point finite -difference approximation provides the following 
relationship: 


ui+l,j + ui-l,j - 2ujj _ 


= (uxx)ij + (uxxxx)jj ^ + o(h^) 


Therefore, with a uniform mesh 


1 

2 



Ui+l,j +'^i-l,j 



( 19 ) 


and overall fourth-order accuracy is achieved. A note of warning should be included 
here. This approximation should not be used with a nonuniform mesh inasmuch as the 
finite -difference approximation introduces a first-order error. Instead, one should 
revert back to the second-order accurate approximation as in equation (17a). 


Examples of Truncation Error Using Burgers' Equation 

Second-order spatial derivative . - With equation (19) approximating Uxx in the 
model Burgers’ equation (7), the system (9), with 9 = 9 \ = 02> takes the form 

AiVf+i^ + = biVj* + + Vf_i) (20) 


If 


Fi = 


At 


2h2 

0 

0 


0 0 

0 0 
0 0 


Gi = 


At 


h^ 

0 

0 


At 


0 

0 


0 

0 


then 


Ai = Ai - 


Bi = Bi + 


Ci = Ci - 
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Di = Di - (1 - 0)Gf 

Ei = (1 - 0)Ff 


It should be noted that equation (19) for hj^ = Constant can be written as 


if * m-l,J - ] + lOMj] + 




“ + 0\h^, 


/ \ ^i-1 i ” ^ ^i+1 1 / a \ 

(Uxx)ij = My + -iiJ J nil * o(h<) 


Note that from equations (18) and (21b) 


(Uxxxx)„ = * oW 


This provides a second-order accurate formula for the fourth derivative. A fourth-order 
finite -difference method developed by H. O. Kreiss is closely related to the present 
spline formulation and is outlined in the appendix. 

If equation (21b) is used for Uxx> then the governing system is still of form given 
by equation (20). The matrices Fj and Gj are now 

0 0 

12 

^i " 0 0 0 (23a) 

0 0 0 


0 ^ 
6 


This paper is concerned primarily with the standard cubic spline approximation as 
given by equations (9). 
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First-order temporal derivative . - The truncation error associated with the time 
discretization in equation (8) or (12) is identical with that found for a typical finite - 
difference formulation. Consider 

ut = f(ui,mi,Mi) 

with 

uP+1 = uP + + (1 - 0)ff|At 

1 1 L -I 

A Taylor series expansion about (n + 0)At leads to 

(ui), = f, + o(!-^ A«t, Al%tt) 

For 0 = 1/2, second -order temporal accuracy is achieved. For all other 0 5 0^1, 
first-order accuracy results. 

For the cases of pure convection {u = 0) or pure diffusion (u = 0) the spline repre- 
sentation for the linear Burgers' equation (7) can be easily transformed into an equivalent 
finite -difference form. 

Diffusion only. - For' u = 0 and i' = Constant, equation (8) with equation (2) can be 
written in the form 

(l - 60|3i+i)un+^ + 2[l + q + 30(/3i+i + u^+l + ^(l - 60^i)un+l 

= [1 + 6(1 - 0)^i+l]uP^l + ai[l + 6(1 - 0)/3[]uP_i + 2[l + Oi - 3(1 - 0)(^i+i + 

(24) 

where (Xi = hj/hi+i and = v At/hf. It can be shown directly from equation (24) or 
by using equations (17a) and (17b) that the truncation error ei is given by 

= utt ^^ ~g^^) At + yuxxxx 22(1 + V) ^ ^ ^t^,(l - CT2)hi+i,hf+i 

For (Ji = 1, the difference equation (24) corresponds to a special case of a more general 
formulation proposed by Saul'yev (ref. 14). Papamichael and Whiteman (ref. 15) recog- 
nized this correspondence in their cubic spline analysis of the one -dimensional heat con- 
duction equation. They considered only the case of a uniform mesh. 

Convection only. - For = 0 and u = Constant, equation (8) with equation (3) can 
be written in the form 
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CT[(l + 30Ci+i)up^j^ + ^2(1 + (Ti) - 30(aiCi+i - CijjuP’*"^ 

= qj^l - 3(1 - 0)ci+ijuf+l 2(1 + CTi) + 3(1 - 0)(cTiCi+i - Cij 


where 


>. - u At 


n 


1 

(25) 


The truncation error is 


®i “ '^ttr^~2/At - uuxxxx 


(1 - 


+ 0 


(At2,h?) 


For a uniform mesh (ctj = 1), equation (25) can be written in the form 


uP+l - uP 
1 1 ~ 

— - + u 


At 


2h 


(ui+i + Uj.i - 2ui)'^+^ - (ui+l + Ui_i - 2uj)" 


uP n - uPl 
+ (1 - 0 ). T +1 


2h 


(6 At )-1 = 0 


(26) 


The fourth-order accuracy is achieved by the effective addition of a difference expres- 
sion representing h^^Uxxt)^- This cancels the u(uxxx)h2 error associated with a cen- 
tral derivative discretization for uux- The largest error terms are now o(h^). 

Fourth-order accuracy with a nonuniform mesh may be possible if a collocation 
procedure is used with a Hermite cubic polynomial approximation. This procedure 
has been analyzed for ordinary differential equations; and fourth-order accuracy can be 
achieved if the collocation points are appropriately located, otherwise only second-order 
accuracy results (see Douglas and Dupont, ref. 11, F5rfe, ref. 13, and Albasiny and 
Hoskins, refs. 16 and 17). 


Complete equation .- For the full Burgers' equation (eq. (8)) a reduction to an equiv- 
alent finite -difference form is possible. Considerable manipulation of the system (9) is 
required, and the final result has not been worked out. The truncation error as obtained 
from equation (11) is 


ei =utt 


1_^20 

e 


r 

j- ^ 


r 

1 

y 

1 + a? 

h?+l - s 

(1 

\ 

12(1 + CT^) 


L 

- 




(1 - crOhihti 


72 


+ 0 


At^, (l - <7j)hi+i,hi+i 
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STABILITY 


General Development for Linearized Burgers’ Equation 

For the linear Burgers’ equation (7) with u and v held constant, interior point 
stability can be assessed with the von Neuman Fourier decomposition of the system (9). 
With 


yn _ yll 0xp Iwx 

or 

vp^e = «” ^ ”1+1 -1 ^ ”l) 

where e = -1, 0, or +1 and I = \/~A, system (9) becomes 

0 


where 


T- 1)0+1 = p.^n 


Ti = 


«0 

^1 


ai 

0 


(27) 


Pq Pi P2 



0 


0 

0 


0 

0 


and the coefficients o!j and pj are defined in equations (9g). Also, 


where 


TTi = hi"^ 1 - exp(-lcpi) + h[+\(l - exp I<Pi+i) 

(28a) 

67T3 = hij^2 + exp(-I(pi)j + hi+i(2 + exp l<Pi+i) 

(28b) 

Tj = 3hj-2 exp [-Icpy) - 1 - 3hj-fi(exp Icp^^i - l) 

(28c) 

T 2 = hi ^[2 + exp(-I<Pij| + hi+\(2 + exp I<Pi+i) 

L J 

(28d) 

= whi. Therefore, 
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uP+1 = GiU^ 

where = Tf^Pj is the amplification matrix. The von Neumann condition necessary 
for the suppression of all error growth requires that the spectral radius 

p(°i) s 1 

The eigenvalues of Gj are found from 
det(Tf^Pi - \il) = 0 
where I is the identity matrix. If 

det Trl = -TTgTg^ao - ^ - «2 = -TrgTgO # 0 

the three roots for \[ are found to be 



For the one -dimensional equation (7), three numerical procedures were considered: 
(1) convection (mi) and diffusion (Mi) explicit, (2) convection explicit, diffusion implicit 
(two steps required for inviscid stability), and (3) diffusion and convection implicit. With 
explicit convection, procedure (1) or (2), both divergence and nondivergence forms of the 
equations were evaluated. 


The stability conditions imposed on these schemes are determined from 



with given by the nonzero value in equation (29). As it is somewhat difficult to eval- 
uate this condition with the expressions (28) for a nonuniform mesh, only the uniform 
mesh stability is discussed here. 


Explicit convection and diffusion . - For a uniform mesh and 0 = 0 in equations (9g), 


1 - 6/3(1 - cos (p)(2 + cos <p)"^ ^ + (3c sin (p)^{2 + cos (p)~^ ^ 1 


where {i = ^ and c = Necessary stability limits are j3 = — , c = (3)”^/^, and 

h^ h 6 

Rc = I = ^ = 2(3)^/^. These results are more restrictive than the limits found for the 
forward time central space explicit finite -difference method, which (from ref. 2) are 
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/3 ^ -i, c = 1, and Rc = 2. In view of this result and the fact that the explicit values for 
mi and Mi must still be determined by the implicit tridiagonal system (2) or (3), this 
explicit representation is not recommended. 

Explicit convection and implicit diffusion . - For = 0 and 02 ~ ^ equa- 
tions (9g), 




1 + 


(3c sin cp)^ 
(2 + cos cp)~^ 


1 + 6/3(1 - cos (p){2 + cos ^)"l] 


-2 


S 1 


This leads to the condition 

c2 s 2^ 
or 



For R(j » 1 this condition is quite restrictive, while for R^. « 1 the stability result 
is quite acceptable. In the inviscid limit Rc the method is unstable as the implicit 
and stabilizing diffusion effect vanishes. This instability can be eliminated if a second 
step is prescribed. This method could then be likened to the Brailovskaya finite - 
difference procedure (ref. 18). The spline technique remains consistent with first-order 
temporal accuracy and second-order spatial accuracy. If the Cheng -Allen viscous cor- 
rection (ref. 19) is made on the Brailovskaya difference procedure, the explicit diffusive 
instability is eliminated; however, the method is no longer consistent in the transient 
unless /3 « 1 (see Rubin and Lin, ref. 20). 

The spline approximation is consistent, and with two steps there is no diffusive 
instability. The two-step procedure is as follows: For step 1, 

u?" = U-* - At(ump + (30a) 

For step 2, 

uP+1 = uP - At(umf + Mf+1) (30b) 

Step 1 is given by system (9) with n+1 — * and 


step 2 is given by system (9) with n - * and 


= 1 a^ = 0 

Pq = 0 Pj = _u At 


Q'2 = At 
P2 = 0 


) 


In addition, the term DVP is added to the right-hand side of equation (9a), where 



0 0 
0 0 
0 0 


( 32 ) 


After Fourier decomposition, step 1 becomes 


Tiuf 


PiU 


n 

i 


p 


where and Pj are defined by equations (9g), (27), and (28), and step 2 becomes 
TiuP"*"^ = P^Uj*' - - v^'j 

or 


where 


Tju 


n+ 

i 


1 


PiU 


n 

i 


Pi = (Pi - D)Ti Ipi + D 

Therefore, 

Pq Pi P2 

0 0 0 

0 0 0 

and 



Pq = 1 + U At yi O' 


Pj = -(u At)2 _1 J 2-1 
2 


P2 = 0 


7^1 

fi = 1 + — j/ At 
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(33a) 
(33b) 


iililiiiiiilbllll il,|yiilllilillilii;.iii 



For a uniform mesh, with subscript i dropped, from equations (28c) and (28d), 

— = -3I(sin (p)\v~^{2 + cos = -Ih“^^ 

"^2 


where I = and 

= 1 + 6i3(l - cos (p){2 + cos 
With Pq’^ 1 replacing Pq,Pj in equation (29) 


|\| 2 = (s2 . c2-i>2)' 



S 1 


If /3 = 0 (inviscid flow), then 


c 5 $ 


-1 ^ 
min 


(2 + cos <p)(3 sin cp) 


-1 


min 


= (3) 


- 1/2 


This result is more restrictive than the c S 1 CFL condition found for the Brailovskaya 
finite -difference method. For 0, the effect of viscosity is to improve the inviscid 
stability limitation; for u 0, the method is unconditionally stable. 

Im plicit convection and diffusion . - For 0 = = 02 0 < 0 = 1 in 

equations (9g), 

1\|2 = _ (1 _ 0)(fi _ i)j2 + (1 _ 0)2g2$2j + 0(j2 _ 1)J2 ^ 02 ^. 2 ^^ < i 

This condition is satisfied and the spline procedure is unconditionally stable if 0 ^ 1/2. 


Development for SADI Procedure 
For the SADI procedure, consider the linear equation 

U^; + UU^ + VUy = ^yy) 

With the spline approximation of equations (14) and a uniform mesh (hij = h and ky = k), 
the amplification factors for the two steps of equations (14) are defined by 


n+ 


2 _ 


GiuU 
1 1] 


„n+l _ 2 

Uij - G2Uij 

where 

Gj. = Tr^Pj. (r = 1, 2) 
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and 


-31 sin (p^ 
6(l - cos ^x) 
-31 sin iPy 
6 (l - cos <Py) 


h(2 + cos <??x) 
0 
0 
0 


h2(2 + cos <Px) 
0 
0 


1 mp m3 


Jij - 


k^2 + cos (Py^ 
0 


k2^2 + cos cp^ 


f} = 

u At 
2 

i\ = . 

< 

1 

II 

COrH 

0 

II 

mj = m^ = 0 


m3 = 

. -V At 
2 

II 

1 

s 

V At 
2 

<2 = 

(| = 0 

()3 _ V At 

2 2 

— 

(1 

-V At 
2 

ml = 

-u At 
2 

m2- 

V At 
2 

"I = “2 ' ° 

— 

<jPx = 

o^xh 

<Py = 

Ct)yk 






The only nonzero eigenvalues of Gj and G 2 are \\ and \2> respectively; that is, 


1 - 3i3y(l - cos ^y)(2 + cos cp^ - 3ICy sin <Py{A 2 cos <p^ 


— 

1 + 3)3^ — cos {2 cos ^ cos ^x) 

1 - 3i3x(l - cos <Px)(2 + cos <Px)”^ - 31cx sin cp-J^A + 2 cos ^^’x)'^ 

^2 = 

1 + 3/3y(l - cos <Py] (2 + cos (p^ "^ + 3ICy sin 9?y^4 + 2 cos (p^ 

where Cy = Cv = ^-r^, and i3„ - From these results it can be 

-*■ h y k h2 ^ k2 
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seen that |^i| |^2| = 1 is always satisfied so that the SADI method is unconditionally 
stable. Of course, boundary effects have not been considered in this interior point sta- 
bility analysis. 


SPLINE CURVE FITTING 


In this section, the accuracy of a cubic spline fit to a given set of data points at pre- 
scribed knots is considered. Error estimates for functional interpolation and for func- 
tional derivatives are reviewed. Exact solutions of the nonlinear Burgers' equation (11) 
are used in a series of numerical experiments to assess the following: (1) mesh require- 

ments and resolution in regions with locally large gradients, i.e., for v « 1 in Burgers' 
equation, (2) the effect of large mesh nonuniformity on overall accuracy, (3) splines under 
tension as a means of smoothing spurious oscillations associated with a cubic spline fit, 
and (4) the accuracy of a spline approximation for the nonlinear convective term (uux)^ 
when this term is obtained from the nondivergence approximation uimi and the diver- 
gence form approximation irij ^where rrq is the spline derivative of the function u^/2, 

i.e., m^ approximates (u^/2)x). This comparison will shed light on the spline solutions 
of Burgers' equation in divergence and nondivergence form. 

Given a set of data points uj at the knots x^ with 
a = xq < xj < xg < < x^^i = b 

if the function u(x), with u(xj^) = uj, belongs to C^ja,!^ then S^(x) in equation (1) 
approximates u(x) to O(h^). Moreover, S^(xj) approximates (aPu/9xP)j to o(h^-P], 
where S]^(xi) = mj and S’^(xi) = (see Ahlberg, Nilson, and Walsh in ref. 10). The 
derivative results have already been inferred by equations (17a) and (17b). If the higher 
order spline approximation, equation (21b) is used to represent (uxx)i> then overall 
fourth-order accuracy for a uniform mesh results. 


Resolution 


A steady -state solution to the Burgers' equation, where = x - Ut, 
ut + (u - U)ujy = VMrjri 


(34) 


with boundary conditions u — 2U as 77 — -oo 

Vv] 


u = U 1 - tan h 


2v 


and u — 0 as 77 — 00 is 


(35) 


Consider the domain -6^77^5. At the knots r) = Vi, u = uj is given by equation (35). 
The spline derivative approximations m| and Mj are found from equations (2) and (3) 
or equation (4). The boundary conditions at 77^ = ±5 are given by 
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TUm^ = 

This represents the spline approximation of equation (34) evaluated at rj = ±5. In sev- 
eral cases, less exact zero slope (mj = 0) or zero curvature = O) conditions were | 

prescribed. These led to negligible changes in the spline curve fit. ! 

The results of a cubic spline approximation as compared with a three -point finite- i 
difference representation, for U=0.5 and selected v values, are presented in tables 1 ' 

to 15. Columns 4 and 8 depict the exact values of found from | 

equation (35). Columns 5 and 9 contain the spline approximations, mj and M^, and | 

columns 6 and 10 depict three -point finite -difference approximations as given by E 

a?Uij_1 - (<J? - llu; - Us 1 " 

(36a) - 


(36b) 

The truncation errors ej for these derivatives are 

h? 

u^ : e^ = + 0(h?^i) (37a) 

2 

^T]r] ■ ®i = (^irTn) ^^ 2 ( 0 -^+ 1 ) ^ 


Wi = 


(" 77 ) i “ 


_ 1 \ 1 


/ ^ 


ai(oi + l)hi+i 


+ 6; 


’’i'^i+1 ■ (^i + l)^i + '^i-]^ 


+ es 


^i(^i + l)hi+l 


For = 0(1), the approximation of equation (36a) is second-order accurate; the approx- 
imation of equation (36b) is first-order accurate unless ctj = 1, in which case second- 
order accuracy results. Therefore, from equations (17a), (17b), (37a), and (37b), S can 
be defined as the ratio of spline truncation error to finite-difference truncation error, so 
that: 




.^ 77777 ) i^i+1 


(" 777 ) i 


{^r]r]T]r}) + 0 

(Ui - l)hf^i,hf^l 


^(^ 777 )^' 

'cf - 1 ) 

(a? + 1 ) 

-1 

- (^ 77777777 )^^ 1+1 + 0 

(a, - 


••77 
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With derivatives of equal order of magnitude and \cri\ S 2, the spline approximation to 
Un should be significantly better than the equivalent finite -difference result. The spline 
fit for Utjtj should be somewhat better for 2 > Uj > 1; however, with a uniform mesh 

{^i = 1 )» 


so that, to lowest order, equal and opposite errors are incurred. This result has already 
been implied in equation (17b) where it is shown that, for a uniform mesh, a spline and 
finite -difference average is fourth-order accurate; this result is depicted in the last 
column of table 4. The exceptional accuracy of the average is apparent. 

In regions of large gradients, with derivatives of increasing magnitude (e.g., bound- 
ary layers or shock waves), the local truncation errors for both spline and finite- 
difference approximations increase. The deterioration is magnified for the spline fit as 
the lowest order truncation error involves higher order derivatives. This difficulty can 
be circumvented with a local decrease in the mesh dimension h^. This mesh reduction 
would be nonuniform so that computer storage is minimized and would be dependent on 
the magnitude of the local gradients. A few additional points properly located can lead to 
a significant improvement in the spline curve fit. The finite -difference formulas are 
also improved, but to a lesser degree. 

For Burgers’ equation (34) when « 1, a region with large gradients, representa- 
tive of a shock wave, develops. The tables for i/ = 1/2 (very weak shock), = 1/8 
(moderate shock), and v = 1/24 (strong shock) show how the spline curve fit varies with 
different placement of mesh points. For the strong shock {u = 1/24) and a uniform mesh 
(h = 0.2), there are few points in the shock region and the derivative approximations are 
poor. With fewer total points but increased density in the shock region, overall accuracy 
is significantly improved. Near the boundaries the derivatives may become smaller than 
the associated truncation errors, and large percentage differences occur. This is par- 
ticularly true with few mesh points. Similar trends can be observed for u = 1/8 or 1/2; 
however, the shock is weaker in these cases, and the agreement is therefore good in 
almost all of the examples presented. In addition to the solutions for a uniform 51 -point 
mesh, only the 15 -point curve fits are depicted. Even for this very coarse grid, the 
spline -derivative approximations are reasonably good for v = 1/8 and much better for 
V = 1/2. With 51 points the agreement is excellent. 

It is significant that for all the cases presented herein, even those in which the 
derivative approximations are poor, the functional valuer between the knots as obtained 
from equation (1) are in excellent agreement with the exact values at the same locations. 
These results reflect the higher order accuracy of the interpolation formula (1). 
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It may be possible to further reduce the required number of knots by using a 
parametric cubic spline curve fit. With such a procedure, step functions can be fit with 
a minimum of knots; this is particularly appealing for regions with very large gradients, 
i.e., strong shock waves. 


Splines Under Tension 

A cubic spline curve fit (eq. (1)) although passing through a prescribed set of data 
points may exhibit spurious oscillations. These oscillations, which are generally much 
less severe than those found with a standard polynomial curve fit, may be suppressed by 
using cubic B-splines (ref. 21), which results in a more complex Galerkin or collocation 
procedure for the solution of differential equations, or by applying tension to the cubic 
spline fit described herein (see refs. 22 and 23). 

In a mechanical sense, tension is used to pull taut the "thread" (curve) passing 
through the data. This results in more accurate interpolation between the knots. The 
cubic spline approximation S^(x) of equation (1) is obtained from a linear distribution 
of the moment S^(x). If a tensile force is added, the spline function in the interval 
[i-l,i], satisfies the following equation: 

^^(x) - (j2s^(x) = (Mj_i - cr2uj_i)(xi - x)hi^ + (Mj - a2uj)(x - Xi_i)hj~^ 

The coefficient o is a tension factor; and for cr = 0, equation (1) is recovered. With 
^a(^ 0 = ^i» = mj, S'^(xi) = Mi, and enforcing the continuity of S)^(xi) so that 

the following results are obtained: 


H,(x) - MiU slnh gfxj - x) Mj sinh g(x - Xj.i) / Mi.A/xi - 
a2 sinh ahj g.2 sinh ahj \ ct 2 jl hs 


Mj\/x - Xj 1 


mi = riMi.i + SjMi + (ui - Ui_i)hfl 


(38a) 

(38b) 


or 


(38c) 

(38d) 


riMi_i + (si + Sj+i)Mi + ri+jMi+i = (ui+j - Ui)h{+\ - (uj - Ui_i)hi'^ 


^i “*’i+l^i+l ■ ®i+l^i + (^i+i - 

mi+l - = (rj+i + Sj+i)(Mi + Mj^i) 


(38e) 



where 


sinh ahs - ah; 
a^hi sinh ahj 

ahi cosh cjhi - sinh ah^ 

= — — 9 

a^hj sinh <?hj^ 


(38f) 

(38g) 


From the relations (38) it can be shown that for splines under tension the coeffi- 
cients of the tridiagonal system (9) become 


ai = ri - ^^j"'^^( 26 iri + 26 i_isi + n)n+l 

bi = Si + Si+i + (- — ( 26i+iri+i + 26iSi+i + 

\hi+lAi+l/ ^ 

- (i^)"^^26iSi + 6i_iri - yi)n+l 
Ci = Ti+i + (- — ( 26 i+lSi+i + 26 iri+i - 

^ \hi+l^i+l/ 


(39a) 


(39b) 

(39c) 


and di is unchanged. For a ^ 0, r^ ^ and si = ^ so that equations (38) reduce to 
equations (1) to (6) and equations (39) reduce to equations (9). For a » 1, there is a 
very large tensile force and equation (38a) becomes 


Sa(x) ^ Ui_i(xi - x)h[^ + Ui(x - Xi_i)hi-^ + 0 (ct-2) 


Therefore Sa(x) is nearly linear between the knots. It has been found (r^f. 23) that, if 
h represents the average mesh spacing, a dimensionless tension factor oh = 1 usually 
works rather well to eliminate oscillatory behavior. 

Examples of curve fitting with tension are given in tables 16 to 24. The tension 
factor CT is generally chosen such that crh = 1. If the mesh resolution is poor (i.e., 
there are insufficient knots in regions of large gradients), tension will smooth the solu- 
tion but accuracy of derivatives is not significantly improved. This plays an important 
role in the solution of differential equations where the functional values at the mesh 
points are unknown and where tension with a poor choice of grid points can lead to a 
deterioration of the solution. This is discussed further in the next two sections. With 
adequate mesh resolution in regions of large gradients, tension does act to reduce oscil- 
latory behavior. 
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Divergence Form 

It is known that for many problems in fluid mechanics accurate numerical solutions 
are possible only when the governing equations are expressed in integral or divergence 
form. For high Reynolds number flows with moderate to strong shock waves and little 
numerical viscosity, numerical solutions obtained with the equations in divergence form 
closely approximate the weak or integral solutions of the differential system (ref. 24). 
This allows for shock capturing by the numerical procedure. Also, for internal flows, 
conservation laws are generally more closely satisfied with the equations in divergence 
form. Bozeman and Dalton (ref. 8) have clearly demonstrated the superiority of the 
divergence form at large Reynolds numbers for the low speed driven cavity problem. In 
regions with moderate gradients or for low Reynolds number flows, there does not appear 
to be any advantage of the more complex divergence -form equations (refs. 1 and 8). 

In order to assess the relative merits and even the necessity for divergence form 
when using splines, curve fits of the nonlinear term in Burgers' equation (34) were 
examined. In divergence form, equation (34) with U = 0.5 becomes 

"* * (^), = 

If mi is the spline derivative of the function and m^ of the function 

u(?7), then the approximation of the nonlinear term is given by m[ in diver- 

gence form and by |ui - ijmi for the nondivergence representation. These expressions 
are presented in tables 25 and 26. Also tabulated are the finite -difference results. 

In the region of large gradients or the shock structure, it is seen that for = 1/8 
the finite -difference approximation in nondivergence form is more accurate than the 
divergence representation. However, for the stronger shock (v = 1/24), the reverse is 
true; and the behavior implies the need for divergence form in high Reynolds number 
finite -difference calculations if thin strong shock waves are to be accurately captured by 
the numerical method (ref. 24). For p = 1/8, the spline approximation provides a sig- 
nificantly better curve fit with nondivergence form of the data; moreover, with divergence 
form, symmetry is no longer maintained. For u = 1/24, this loss of symmetry is still 
evident, but neither of the spline curve fits is particularly good as long as the grid den- 
sity in the shock region is too low; however, it will be observed that for the spline form- 
ulation the solutions in nondivergence form are generally as good as or better than those 
obtained with divergence form. These results agree with the conclusion of Douglas and 
Dupont (ref. 11) that, unlike Galerkin or finite -difference methods, there appears to be no 
advantage to divergence form for collocation procedures. 
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BURGERS' EQUATION 


Numerical Procedure 

The spline calculation procedure was first applied to the one -dimensional nonlinear 
Burgers' equation (34) with boundary conditions u 2U as 77 "°° u 0 as 

77 — <». Initial conditions were specified such that 

^0 (v7 > 0) 

u(t7,0) = ( U (77 = 0) 

2U (77 < 0) 

For all of the solutions presented here, U = 0. 5, v = Constant, and the boundary condi- 
tions are prescribed at the finite locations 77 = ±5. 

The tridiagonal system (10) is used to determine MP+1; then uP+1 and mP+1 
are obtained from equations (5), ( 6 ), and ( 8 ). The predicted stability restrictions were 
all confirmed by the calculations. Therefore, the solutions to be discussed here were 
obtained either with the implicit (9 = 1) procedure or the two-step method. The former 
is unconditionally stable; the latter has the restriction c ? 

Implicit method. - For the implicit method, the nonlinear coefficient 26j = u[ At 
= (ui - 0 . 5 ) At in system (10) was linearized by evaluation at the latest time t = n At. 
Since the steady-state solutions were obtained in a minimum of computer time, the more 
accurate quasi -linearization procedure of equations (11) to (13) was unnecessary. The 
convective derivatives as found from equations (5), (6), and (8) are given by 


mP-^1 = - 


6^ + 26j_i 

3A^ 


hiAi 




\3 


3A; 


hiAi 


Uj - u 


i-1 


Mi 


The tridiagonal system (9), although unconditionally stable, is only conditionally diagonal 
dominant. For a uniform mesh, diagonal dominance is assured if 


Rc 


_ c _ uh < o 


For Rp > 2, diagonal dominance requires 



or for Rq — 00 , 



A similar result is found with finite differences (see Hirsh and Rudy, ref. 25). 
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A correction similar to that proposed by Khosla and Rubin (ref. 26), which provides 
diagonal dominance of an implicit finite -difference method, can be formulated for this 
spline procedure. In view of the fact that the tridiagonal system (9) was inverted for c 
as large as 600, this correction was unnecessary for the present calculations. The final 
solutions were invariant for ^ ^ 600; this corresponds to a maximum of 431 time 

steps to convergence and to a minimum of 22 for the conditions: 31 points, v = 1/24, 
unequal spacing, = 1.5, and a = 0. 

Two-step method . - The two-step method is outlined in equations (30) to (32). The 
tridiagonal system (10) is obtained for each step, with 0 = 0 for terms and 0 = 1 
for 7 | terms. Since the convective terms are evaluated explicitly in each step, solu- 
tions were obtained with both divergence and nondivergence forms of the convective deriv- 
ative. In nondivergence form. 


(26imi)’^ = (ui - 0.5)^mp At 

with mP obtained from equation (5) or (6). In divergence form, 
(26imi)" = AP At 


where mP is obtained from equation (3), with m^ — in^ and uj — 


Ui 


The boundary conditions for both the implicit and two-step procedures, applied at 

(for V * 1/2, the exact values 0.9933 and 0.0067 were 
The procedure outlined previously leads to a linear 
Mjvf and respectively. In several cases, the 

Mj = 0 were tested; these did not cause any significant 


m = ±5, were uq = 1, = 0 

prescribed), and T0.5m = 
relation for Mq and Mj and 
bovindary conditions mj = 0 or 
variations in the solutions. 


Numerical solutions were also obtained for splines under tension as described by 
the systems (38) and (39). The procedures for the implicit and two-step methods are 
identical with those of the preceding discussion. 


Discussion of Numerical Results 

The steady-state results of calculations for Burgers' equation are given in tables 1 
to 24 and for the strong shock {v = 1/24) in figures 1 to 17. The two-step solutions are for 
divergence form if unspecified. In the tables, columns 3, 7, and 11 depict the calculated 
values of ui, mj, and M^, respectively, as obtained with the implicit nondivergence - 
form spline procedure. The two-step results in nondivergence form were almost identi- 
cal with the implicit solutions. This can be seen from several of the figures where both 
solutions are depicted. 
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In table 9, column 12, the calculated values of U[ with the two-step divergence - 
form procedure are presented. These are considerably less accurate than the implicit 
solutions and follow the pattern of the curve fits in tables 25 and 26. This loss of accu- 
racy with divergence form is found for a nonuniform mesh as well and leads to a conclu- 
sion that is contrary to that of finite -difference procedures. For spline calculations, the 
nondivergence form appears to be preferable providing there is adequate grid resolution 
in the shock structure. This result may be significant for any shock-capturing study. 

The solutions for u[, as seen from the tables and figures, are all quite acceptable. 
This is true even when the agreement for m| and Mj, with the analytic derivatives 
obtained form equation (35), is not very good. As the number of mesh points in the shock 
structure is increased, a noticeable improvement in m[ and Mi can be discerned. 

This behavior follows the trend found for spline curve fitting; i.e., accurate derivative 
representation in regions of large gradients is possible only with adequate mesh resolu- 
tions in these regions. In this respect, a nonuniform mesh requiring fewer total grid 
points is desirable. As few as 15 to 19 points leads to solutions for uj that are quite 
good; moreover, as seen with = 1.6 and 19 total points, in the shock ~ 0.044, 

while near the boundary = 1.905. This represents a significant mesh variation, 

which is highly desirable fdr shock and boundary -layer problems. 

Another interesting feature of the spline solution is shown in figure 1, where a non- 
divergence, central derivative, finite-difference calculation is also' depicted. Nondiver- 
gence finite -difference solutions for a nonuniform mesh are shown in the last column of 
tables 5 and 6. The difference formulas are given by equations (36) and (37). Since 
(^c)max ^ finite -difference solution exhibits an oscillation typically 

found for Rc > 2. This oscillation does not occur with the spline result and may be indi- 
cative of the fourth -order accuracy of the convection derivative mj. In addition, the 
interpolation formula (eq. (1)) provides intermediate u{r]) values that agree very well 
with the exact solution (eq. (35)) so that accurate grid realinement is possible. 

In several runs, with a nonuniform mesh and few mesh points, oscillations did 
appear in the final solution. In most cases, tension (with oii = 1) removed or minimized 
this spurious behavior; h is the average mesh width (lO/N+2). The primary effects of 
tension can be summarized as follows: First, for the 15-point calculations with boundary 

conditions at 77 = ±5, hj near the boundary is extremely large, e.g., (hi)j„ax = 2.261 
with CTj = 1.8, so that stable solutions are obtained only when tension is included in the 
spline formulation. Second, for 15 points with boundaries at 77 = ±3 or 19 points with 
boundaries at 77 = ±5 stable solutions are obtained, but oscillations appear near the 
boundaries. Tension has the effect of minimizing or eliminating the oscillations with no 
apparent loss of accuracy. Third, if the spline derivatives are inaccurate, as with a 
coarse mesh in the shock structure, tension does not improve the accuracy but appears 
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to have a negative effect. This is apparently due to the low mesh density in a region of 
large gradients; a similar effect occurred with the spline fitting under tension as previ- 
ously discussed. Fourth, it does appear that tension can be used to smooth oscillations 
arising from large changes in mesh width when local singular regions form and increased 
mesh resolution is required. 


DIFFUSION EQUATION 


The two-dimensional diffusion equation 
"t = g(uyy + Uzz) 


where u = u(t,y,z), with the initial condition u(0,y,z) = 0 and with the boundary condi- 
tions u(t>0,0,z = 0) = 1, u(t>0,y=0,0) = 1, and u(t,y,z) -*■ 0 as y,z — °o has the exact 
solution 

u = 1 - erf Y erf Z 


where ^ ~ 


This solution describes the impulsive motion of 


a right-angled corner formed by two infinite flat plates and has been used by Sowerby 
(ref. 27) to infer the steady flow along the corner with leading edge at t = 0, i.e,, 
Rayleigh's problem for a corner. This problem was used to test the accuracy of the 
SADI procedure previously outlined. The two-step procedure is given by 



u 


+ 4i L; 




and 


1 

-r' = "li ' 



+ 



where and Pjj are the spline approximations to (uyyjjj and (uzz)jj, respec- 

tively. The boundary conditions are simply uy = 1 on the walls and uij — 0 as 
y,z — oo. In addition, from the governing equation, Ljj = 0 on y = 0 and z > 0 and 
Pjj = 0 on z = 0 and y > 0. 

Some results of this SADI calculation for R = 1000 are given in table 27 and fig- 
ure 18. A nonuniform grid was specified in order to accurately describe the boundary- 
layer behavior near the walls with a minimum of mesh points. The agreement with the 
exact solution is reasonably good so that the validity of the SADI procedure is confirmed 
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III \m. 



INCOMPRESSIBLE FLOW IN A CAVITY 


Numerical Procedure 

As a final test problem the incompressible flow in a driven cavity was considered. 
This problem has been studied by numerous investigators, and recently Bozeman and 
Dalton (ref. 8) have reviewed the literature and presented some definitive results and 
conclusions. The governing equations in terms of a vorticity stream -function system 

are 


*^xx + ^yy ^ 


(40a) 


h + ^^x + '^^y " r(^xx + 


^yy) 


(40b) 


where 4/ is the stream function, C is the vorticity, and u = i//y and v = are 

the velocities in the x- and y-direction, respectively. The boundary conditions and 
geometry are shown in figure 19. For all the calculations the initial conditions are 

i//(x,y) = 0 and ^(x,y) = 0. 

Solutions are obtained by an iterative SADI procedure. The SADI system repre- 
senting equation (21) is given in two steps for both stream function and vorticity. 

Stream function. - For step 1, 


n+l,s+4 , „ . 

^i. 2.^n.l,s,^r 




n+l,s+i 


(“«) 




- C 


1] 


(41a) 


For step 2, 


^n+l,s+l = 

^1] 1] 


n+l,s+4 


At 


Li^' 

1] 


n+l,s+^ 




^^n+l,s+l _ 




1] 


(41b) 


The physical time t equals n At; At is the time increment at each step n; At is 
a fictitious time step; and t = s At. Solutions for equation (40a) are obtained as the 

steady-state limit (t -«)) of equations (41); Ly and Mij are the spline approxima- 
tions to ^ and respectively. (The superscript 4y implies that A = i//.) 

First derivatives m and = -v are represented by and my, respectively. 


Vorticity. - For step 1, 


n+4 

V 2 _ >»n At 
2 




"4 


n+i 


9 


(42a) 
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For step 2, 


i»n+l 








+ R-1 



+ R-1 



(42b) 


The bar over the spline derivatives denotes an average of n and n+1 values; the 
? superscript denotes ^ spline derivatives. 


The iterative procedure is as follows: 

(1) Given i//jt and either as initial conditions or at time n At, all the 
spline derivatives are determined from equations (14) to (16) and (4) to (6). On the ver- 
tical surfaces, M^j = ?ij; on the horizontal boundaries, Ljj = 

(2) The vorticity is obtained with the SADI technique as outlined in equa- 

tions (42) and (9). At the boundaries, is found from an expression similar to equa- 
tion (5) or (6). At the upper moving wall (w), with 


ip _ _ i^iw^iw •^iw^i^w-l (’^iw - ’^i,w-l) 

3 6 ki 


'^iw 


and with = 0 and the vorticity becomes 




r. = ^ _ 
^iw 


^i,w-l 


k; 


iw 


k? 


Similar relations can be derived for the three stationary walls. Boundary values for 
M^j and are obtained from equations (42) evaluated at the surface. For the moving 

t» n+ — 

wall, equations (5) and (6) are used to eliminate my. In addition, Ciw^ is evaluated 
with the three -point formula 


n+-i 


Cl 


2 = 


35";* + 6C"w - 


aw 


IW 


0(At2) 


(3) The vorticity is used in equations (41) and the SADI procedure is applied 

over the fictitious time t until a converged solution, to any specified tolerance, is 
obtained. 


(4a) If only the steady-state solution is required, the calculation proceeds to the 
next time step (n+2) by returning to step (2) with n - n+1. The spline derivatives for \p 
have already been determined in step (3). 
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(4b) If an accurate transient is required, the calculation proceeds to step (2) with 
and all spline derivatives of ip replaced by averages over the n and n+1 time 
steps. Then is recalculated, and this process continues until convergence. 

Although accurate transient solutions have been obtained in a number of cases, only 
the steady-state results are presented here. The time step At was generally chosen 

such that At = (h;4,kii) . Larger values were used in many cases, but a careful study 

\ J J/ min j. x*n v,. 

of optimal time integration by discrete or semidiscrete procedures must still be consi- 
dered. Primary interest at this time was concerned with the applicability of the SADI 
procedure, as well as the accuracy, ease of handling boundary conditions, and other gen- 
eral characteristics of the spline approximation. 

Calculations are presented for a square cavity with R = 10 and 100 and for a rec- 
tangular cavity with b/a = 2 and R = 100. Comparisons are given with finite -difference 
calculations in both divergence and nondivergence form. Central differences are used 
throughout. The vorticity equation is solved with an ADI procedure, and the solution for 
i// is obtained by a direct Poisson solver or by successive overrelaxation. 

Discussion of Numerical Results 

The results are presented in tables 28 to 34 and figures 20 to 22. For all cases, 
the values of i^'niax vorticity ? at the midpoint of the moving wall are 

depicted. For these values, comparisons between the spline and finite -difference solu- 
tions are possible even when the grid alinements differ. In addition, the distributions of 

xp/ and ? for the spline solutions, and in several cases for the finite -difference solu- 
tions, are presented. The figures depict the horizontal velocity component, uj^j or 
f jj, along a vertical line passing through the vortex center. 

The results for R = 10 are given in tables 28 and 29 and figure 20. For this low 
Reynolds number the spline and finite -difference solutions in either divergence or nondi- 
vergence form are quite similar. For R = 100 the large disparity between the diver- 
gence and nondivergence finite -difference solutions, first noted by Bozeman and Dalton 

(ref. 8), is apparent. The values of and ?^all ^ 

variety of grids. Also included is a limiting solution obtained by Richardson extrapola- 
tion (ref. 12) from the two or three calculated values of each procedure. It is evident 
that the divergence finite -difference solution is more accurate than the nondivergence 
result; however, the spline solution, which is obtained in nondivergence form, appears to 
be even more accurate than the divergence -form finite -difference result. For example, 
the value of <^'jjiax obtained from the spline calculation with 15 points (this denotes 
a 15 ^ 15 node mesh with h = k = 1/14) is about 1 percent higher than the extrapolated 


35 


value; the 17-point divergence -form finite -difference result is about 4 percent lower 
than its extrapolated value. The nondivergence 15-point finite -difference result is low by 
about 12 percent. These results again seem to reflect the higher order accuracy of con- 
vection terms in the spline procedure; in the vortex core region, the flow is inviscid 
dominated. However, near the moving wall, where diffusion is important, the vorticity 
results appear to show a similar trend; the spline values are always somewhat more 
accuratei than the divergence finite -difference solutions. 

Another interesting result is shown in table 30(c). The velocity at the first grid 
point away from the upper left boundary is depicted. With 15 points, the spline result is 
of an opposite sign to that obtained with the nondivergence finite -difference method. With 
a finer grid, the finite-difference solution changes sign so that once again the spline pro- 
cedure prevails. Unfortunately, the more accurate divergence -form finite -difference 
solutions were obtained with a slightly different grid so that a direct comparison is not 
possible. However, a change in sign with mesh reduction is observed for the velocity in 
the corner region. An interpolation procedure is used to estimate the value at the 
desired location. These results are also given in table 30(c). The extrapolated limit 
closely approximates the solution obtained with splines. The velocity profiles through 
the vortex center are shown in figure 21. These values are also tabulated in table 31. 

The agreement is quite good. 

A spline solution for R = 100 was also obtained with a 19 -point nonuniform mesh. 
In the central region of the cavity hy = kjj = 1/14 as with the 15-point mesh; however, 
near the boundaries there is some grid realinement to Increase the mesh density in the 
surface boundary layers (see table 32(g)). The increased accuracy near the boundaries, 
where diffusion is most important, leads to a solution that appears to be almost as accu- 
rate, throughout the entire flow domain, as the 29-point results. The improved accuracy 
of this 19-point solution is seen in tables 30(a) and 30(b) where </^niax ?wall 

indicated. These results imply the considerable advantages of the spline procedure with 
a nonuniform grid in regions of large gradients. In this manner, the accuracy of the 
second-order diffusion terms is enhanced in domains where these effects are significant. 
In inviscid regions the fourth -order accurate convection terms are dominant, and mesh 
reduction is not as important. The improved resolution of the corner vortices is seen in 
the distributions of tables 32(g) and 32(h); the comparisons with the 65-point diver- 

gence finite -difference solutions are reasonably good. 

For R = 100, spline solutions were also obtained for a 2 x 1 rectangular cavity with 
a 29 X 15 point uniform mesh; the results are presented in tables 33 and 34 and figure 22. 

A double vortex is observed. The flow properties are in qualitative agreement with the 
divergence -form finite -difference solutions obtained with a 33 x 17 uniform mesh. 


36 


CONCLUDING REMARKS 


The use of a cubic spline approximation for the evaluation of spatial gradients pro- 
vides a highly efficient and accurate procedure for numerical calculations with a uniform 
or nonuniform mesh. It has been shown that: (1) Second-order spatial accuracy is 

achieved, even with an arbitrary nonuniform mesh, for equations of the Navier -Stokes 
type; (2) For inviscid regions, with a nonuniform mesh, third-order accuracy results; 

(3) For the Navier-Stokes equations and uniform mesh, the interior point truncation error 
is fourth order with a combined spline finite -difference scheme; (4) Derivative boundary 
conditions can be treated easily and accurately so that spatial finite -difference discreti- 
zation is unnecessary; (5) There appears to be no particular advantage gained with the 
divergence form of the equations; (6) Accurate interpolation is possible if grid realine - 
ment becomes desirable; (7) Evaluation of quadratures which are generally not of a tri- 
diagonal form, as in finite -element or other Galerkin procedures, is unnecessary. 

With a finite -difference discretization for the time -like integration, it has addition- 
ally been shown that: (1) The system of algebraic equations resulting from the spline 

formulation is block tridiagonal, and therefore inversion for implicit time discretization 
is accomplished with an efficient algorithm. Moreover, appropriate substitutions can 
reduce the vector system to a scalar one, thereby eliminating the necessity for any 
matrix inversions. (2) Explicit, implicit, and mixed time integrations have been consid- 
ered. The interior point stability conditions for explicit procedures are slightly more 
restrictive than those found with equivalent finite -difference techniques. Implicit methods 
are unconditionally stable. 

Solutions have been obtained for the one -dimensional nonlinear Burgers' equation, 
and in two dimensions for the diffusion equation and the vorticity-stream function system 
depicting the incompressible viscous flow in a driven cavity. Oscillations typically found 
with second-order accurate finite -difference methods when the cell Reynolds number 
exceeds 2 did not occur with the spline solutions; and this probably reflects the higher 
order accuracy in the convective term. Accurate solutions for Burgers' equation are 
obtained even with a highly nonuniform mesh if adequate mesh resolution is specified in 
the region of largest gradients. 

For two-dimensional flows the SADI procedure appears to work quite well for both 
the diffusion equation and driven cavity problem. Comparisons with the analytic solution 
available for the former are excellent, and with finite -difference calculations for the lat- 
ter are quite reasonable. The spline solutions for the cavity obtained with the non- 
divergence form of the equations are somewhat better than the divergence form finite - 
difference solutions and considerably better than the nondivergence form finite -difference 
results. Once again the higher -order accuracy of the convection operator may account 
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for the improvement with the spline formulation. The vorticity boundary condition has 
been treated directly and without the need of any finite -difference discretization at the 
boundaries. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., January 30, 1975. 
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APPENDIX 


KREISS FOURTH-ORDER FINITE-DIFFERENCE METHOD 

For a uniform mesh, H. O. Kreiss has proposed a fourth-order method that is very 
similar to the basic spline procedure presented here, see Orszag and Israeli (ref. 28). 
For Burgers' equation (7), Kreiss' method reduces to the system of equations ( 8 ), (2), 
and (3) except that the coefficients h/ 6 , 2h/3, and h /6 for Mi_j, Mj, and 
respectively, in equation (2) become h/12, 5h/6, and h/12, respectively. No longer 

is Mj a spline approximation but a finite -difference approximation such that 

Mi = (uxx)i + O(h^) 

The system (9) describes Kreiss' method, with hi /6 — h/12 in Ai, 

and i in Cj. All other entries in equation (9c) are unchanged. 

O 1a 

The stability of this procedure can be assessed directly from equation (29); and 
Pj are given in equations (9g); tj, T 2 , and ttj are given by equations (28). Due to the 
change of coefficients in equation ( 2 ), 

6173 = h(5 + cos (p) 

instead of the spline value 

6773 = h(4 + 2 cos 0) 

The stability condition |\i| ? 1, with \i given by the nonzero value in equation (29), 
leads to the following results; (1) The implicit procedure (0 = 1) is unconditionally 
stable; and ( 2 ) the explicit procedure (0 = 0 ) has a stability condition 

/l . 12/! . |'3csin_^/ , J 

\ 5 + cos 0/ \2 + cos 0/ 

Therefore, necessary stability restrictions are jS = 4, c = — , and Rg S 

o vG ^ 
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TABLE 2 .. IMPLICIT SPLINE SOLUTION TO BURGERS' EQUATION FOR >-=1/2, 0 = 0, q = 1.4, AND 15 POINTS 





TABLE 4.- IMPUCTT SPLINE SOLUTION TO BURGE2RS' EQUATION FOR u = 1/8, ct = 0, AND 51 EQUALLY SPACE PCXNTS 














TABLE 5.- IMPLICIT SPLINE SOLUTION TO BURGERS’ EQUATION FOR u- 1/8, a = 0, aj = 1.2, AND 15 POINTS 
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TABLE 7.- IMPLICIT SPLINE SOLUTION TO BURGERS’ EQUATION FOR v = 1/8, = 0, (tj * 1.6, AND 15 POINTS 
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TABLE 9,- IMPLICIT SPLINE SOLUTION TO BURGERS’ EQUATION FOR u = 1/24, a = 0, AND 51 EQUALLY SPACED POINTS 
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TABLE 14.- IMPLICIT SPLINE SOLUTION TO BURGERS’ EQUATION FOR i/ * 1/24, tr = 0, C7^ = 1.4, AND 19 POINTS 



Spline 

calculated 

M 

eocacaeqc«*Hi -4 

0000000 0 

yyyyyyy y 

eQ*-*poo>ooo»«eooocsi 

«' 4 *pao<ocaTi*oco« 

eoo'^N'-'OOeotrtTr 

oe >-< *-< ci V C-’ TT csi 

1 1 t t T-l 

Finite - 
difference 
curve fit 
U777? 

1 CO 

I •-< os 0 Tf « 

1 ) ( 1 1 t ) 

1 ® 2 2 2 2 S 

{ y y y y y y 

j ifsirsost-ooO'^Os 

j ocamrfosoost- 

1 w CO 'st* os ^ os ifs »-• 0 

1 1 1 1 t 1 > 1 w 

Spline 
curve fit 
M 

Tp CO 

caift-o*eocoC 3 T-t t-i 

1 1 1 1 1 f t 1 

0000000 0 

xyyyyyy y 

Wt-OOOO'^COOS—lC- 
oo«oo 5 cooee^i>ooco 
c« 3 eoirttft»HQocomos»-i 
1-4 r-' co’ 1-4 r-’ eo 1-^ in id i -4 

II t 1 1 I --I 1 

a " 

. ^ 

« ^ ^ r- ^ ^ 

II 1 1 1 1 1 

0000000 
y y y y y y y 

i-llD‘^^C<li-<{»C-’^eM 
oesit-i-f-^raT-^ooe 
pfliDcococomowt- 
w 00 CO co’ ca N t-' 0 ed 0 

1 1 1 1 1 1 1 1 l-H 

Spline 

calculated 

!: 1 

■"U'-O* cocoeoeOMi-t 
II 1 1 I 1 1 1 

00000000 

yyyyyyyy 

<^OON®OSPOC 40 '^ 
TpTpoococooincO'— <LO 
OS«Oi-iOO' 5 *if’coUOCJO> 

CD ® 1-4 1-4 ed ® ed to cd cd 

1 1 I t 1 1 1 

Finite- 
difference 
curve fit 

U 77 

« eo 

I 1-1 os to iir N 1-t 1-f 

II 1 I 11 1 1 

1 0000000 

I y y y y X y y 

' rHCii#<irt®ot-eo® 

if<i/scoi-ii-ieotomc- 

j N to 0 Tf M 0 to 

1 .-4 .-4 ,-4 ed 1-4 cd ed ed 

1 1 I I I 1 1 1 1 t 

Spline 
curve fit 
m 

uomTt»ij*eococ^i-f 

00000000 

yyyyyyyy 

^tOc^eMCMTPOt-®® 

■^OONlOMOSOeO'iJi 

COtDi-*tO’-^OU^O'-<OS 

»4 ed -4 ed id ed to’ c 4 ed 

t> 

»o ® ca 

caw i-icoiocscaw 
00000000 

yyyyyyyy 
woeoowootoioc-o 
mocwcoiooOcrjwo 
owoor-oswosow O 
w r 4 ci ci -4 cd id id ca ed 
1 1 1 1 1 1 1 1 I 1 

T 3 

c is 

n 

Oc^weoi-cffgcatDijiO 

OOOOOOOtOt-O 

0 0 0 0 0 0 0 os t- ® 

W 1— 1 W »H W W 1—1 

0 

rt 3 
& 

OOoOOOifiOSWO 

ooooooos^t^o 
0000^ 0 0 OS OS t- in 

R- 

0000 
y y y y 

oioeot-wioocaco 
oosc 3 inwocae»sw 
o^-^eOwdif^o 
id cd <d w 1-4 t 4 'id ed 1-4 0 

1 1 I 1 1 I 1 1 I 


52 




TABLE 16.- IMPLICIT SPLINE SOLUTION TO BURGERS’ EQUATION FOR = 1/8, a = 2.0, q = 1.4, AND 15 POINTS 
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TABLE 22 ,- IMPUCTT SPLINE SOLUTION TO BURGERS’ EQUATION FOR ^ = 1 / 24 , o = 5 . 0 , = L 4 , AND 19 POINTS 


Spline 

calculated 

M 

0 

in LO ^ rr « CO N <1* 

0000000 0 

XXXXXXX X 

mocoOtoM^iC^O 
00 « *« ® 0 r- 0 0 VI 0 

M vS t-J CO •-< «n « eo vi 

Finite - 
difference 
curve fit 

'ijJTJ 

1 CO ' ^ 

1 >-( 0 <0 O' csi rH 

t 000000 

1 X X X X X X 

• mmoi-^WO'^o 

' — <tt),-(moo-a*t-co 

1 »o « 0> 0 0 

1 ,-ieO'^O'^OUi^O 

Spline 
curve fit 
M 

eo 

com^TfeoWri >-* 

0000000 0 

XXXXXXX X 

<— <CD^O'u 3 cocoOC'*-^ 
coVicoo- 4 '’"*oo<fi*rJ 
M^rcot— iAovioco^ 

wr-Ii-Iijivicociioinr-i 

Exact 

N ^ c- ^ 

0000000 
XXXXXXX 

r-<c 3 tr)C'j'-<CDr-^cc 

CflO«OCOCOV 50 Ni> 

»-(' 00 CO CO cvi Cfl cc CO 0 

Saline 

calculated 1 
m 

50 cotomto^co^ 
00000000 
XXXXXXXX 

coc-'-^OOS'-'co-'J'r-oj 

oeocoNC'^^-o^^^ 
eo t- 0 « 0 c-^ 0 00 

c<ici-^cioco-^t^C 3 Cji 

Finite - 
difference 
curve fit 
U77 

1 CO 

1 04 m ^ 

1 0 0 0 0 0 0 0 

; X X X X X X X 

' r-tC^^mODOF-COOO 

' Ti>inco>-<’-'co<naDr- 

' oai>-r-tDO^CflOtD 

Spline 
curve fit 
m 

t- (O V 5 ^ Tf CO 

00000000 

XXXXXXXX 

CCOT-HTT 04 -- 'COOCTfC— 

oe lO Neoco N 0 

CO pa ca •'ji 50 ^ *? *?*?'*? 


tn » N 

pa r-iwoomcoPi*-* 
00000000 

XXXXXXXX 

r-«coeoo«^® 5 DV 5 r-o 

in®r^tDV>ooco<-«0 

o«Hoot-a»*^04®'“50 

•Hc-paca'-ipatniopieo 

Spline 

calculated 

g8S8SSgS2g 

00000000 t-^V 5 

Exact 

u 

OO 00 QOV 5 O -^0 

Ooooooo-^r-o 

oooooooot^vi 


0000 
X X X X 

05 DP 4P"’“'V*OCa® 
oopaio’-'opacow 
0 ^ 'J* *0 ’-J w 0 

inco'eJ'-«'-'P-^N'“< 0 

. ’ ' ’ ' 




57 


TABLE 25.- COMPARISON OF SPLINE AND HNITE- DIFFERENCE CURVE HTS OF THE EXACT SOLUTION 
TO BURGERS’ EQUATION FOR i' » 1/8, a => 0, AND 51 EQUALLY SPACED POINTS 


V 

Exact 

u 

Spline 
curve fit 
(u-l)m 

Spline 
curve fit 
m 

Exact 

Finite -difference 
curve fit 

("-jK 

Finite-difference 
curve fit 
d u2 1 

-5.0 

-4.8 

0,9999999980 

.9999999950 

0 

-8.91817 X 10-6 

0 

-1.07708 X 10-8 

-4.12231 X 10-9 
-9.17436 X 10-9 



-1.01850 X 10-8 

-1.01850 X 10-8 

-4.6 

.9999999900 

-2.04316 X 10-8 

-2.59096 X 10-8 

-2.04179 X 10-8 

-2.26687 X 10-8 

-2.26688 X 10-8 

-4.4 

.9999999770 

-4.53333 X 10-8 

-5.76708 X 10-8 

-4.544 09 X 10-8 

-5.04462 X 10-8 

-5.04463 X 10-8 

-4.2 

.9999999490 

-1.00921 X 10-7 

-1.28429 X 10-7 

-1.01131 X 10-7 

-1,12267 X 10-7 

-1.12268 X 10-7 

-4.0 

.9999998870 

-2.24592 X 10-7 

-2.85815 X 10-7 

-2.25070 X 10-7 

-2.49856 X 10-7 

-2.49856 X 10-7 

-3.8 

.9999997500 

-4.998 54 X 10-7 

-6.36108 X 10-7 

-5.00903 X 10-7 

-5.56073 X 10-7 

-5.56074 X 10-7 

-3.6 

.9999994430 

-1.11242 X 10-6 

-1.41567 X 10-8 

-1.11478 X 10-8 

-1.23755 X 10-8 

-1.23755 X 10-5 

-3.4 

.9999987600 

-2.47574 X 10-5 

-3.15061 X 10-8 

-2.48098 X 10-8 

-2.75421 X 10-8 

-2.75420 X 10-5 

-3.2 

.9999972390 

-5.50985 X 10-6 

-7.01175 X 10-8 

-5.52148 X 10-6 

-6.12957 X 10-8 

-6.12956 X 10-6 

-3.0 

.9999938560 

-1.22623 X 10-5 

-1.56046 X 10-5 

-1.22881 X 10-5 

-1.36414 X 10-5 

-1.36413 X 10-5 

-2.8 

.9999863260 

-2.72893 X 10-5 

-3.47272 X 10-5 

-2.73469 X 10“5 

-3.03584 X 10-5 

-3.03581 X 10-5 

-2.6 

.9999695680 

-6.07295 X 10-5 

-7.72793 X 10-5 

-6.08576 X 10-5 

-6.75586 X 10-5 

-6.75572 X 10-5 

-2.4 

.9999322760 

-1.35136 X 10-4 

-1.71951 X 10-4 

-1.35421 X 10-4 

-1.50328 X 10-4 

-1.50321 X 10-4 

-2.2 

.9998492900 

-3.00652 X 10-4 

-3.82503 X 10-4 

-3.01284 X 10-4 

-3.34432 X 10-4 

-3.34398 X 10-4 

-2.0 

,9996646500 

-6.68627 X 10-4 

-8.50373 X 10-4 

-6.70026 X 10-4 

-7.43649 X 10-4 

-7.43481 X 10-4 

-1.8 

.9992539710 

-1.48564 X 10-3 

-1.88809 X 10-3 
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5.57393 X 10-'^ 

1.11243 X 10-8 

1.15016 X 10-8 
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2.49856 X 10-7 

4.2 

5.05660 X 10-6 

1.00922 X 10-7 

1.04344 X 10-7 

1.01131 X 10-7 

1.12267 X 10-7 

1,12267 X 10-7 

4.4 

2.27200 X 10-6 

4.53 286 X 10-8 

4.68685 X 10-8 

4.54409 X 10-8 

5.04462 X 10-5 

5.04462 X 10-8 

4.6 

1.02090x 10-8 

2.04423 X 10-8 

2.11259 X 10-8 

2.04179 X 10-8 

2.26687 X 10-8 

2.26687 X 10-8 

4.8 

4.58500 x 10-9 

8.91511 X 10-9 

1 9.07434 X 10-9 

9.17436 X 10-9 

1.01850 X 10-8 

1.01850 X 10-8 

5,0 

j 2.06100 X 10-6 

L» 

L“ . 

4.12231 x 10-9 




TABLE 26.- COMPARISON OF SPLINE AND FINITE -DIFFERENCE CURVE FITS OF THE EXACT SOLUTION 
TO BURGERS' EQUATION FOR v = 1/24, a = 0, AND 51 EQUALLY SPACED POINTS 


rj 

Exact 

u 

Spline 
curve fit 
(u-l)m 

Spline 
curve fit 
m 

Exact 

(“-ih 

Finite -difference 
curve fit 
(u-^u. 

Finite-difference 
curve fit 

A U? . 1 u„ 
drj 2 2 ^ 

-5.0 

1,0 

0 

0 

-5.25391 X 10-28 



-4.8 

1.0 

1.27547 X 10-14 

-1.12666 X 10-14 

-5.79147 X 10-25 

0 

0 

-4.6 

1.0 

-4.46447 X 10-14 

3,94928 X 10-14 

-6.38404 X 10-24 

0 

0 

-4.4 

1.0 

1,65835 X 10-13 

-1.46687 X 10-13 

-7.03724 X 10-23 

0 

0 

-4,2 

1.0 

-6.18734 X 10-13 

5.47288 X 10-13 

-7.75728 X 10-22 

0 

0 

-4.0 

1,0 

2.30925 X 10-12 

-2,04260 X 10-12 

-8.55098 X 10-21 

0. 

0 

-3.8 

1.0 

-8.61880 X 10-12 

7,62358 X 10-12 

-9.42590 X 10-20 

0 

0 

-3,6 

1.0 

3.21680 X 10-11 

-2,84535 X 10-11 

-1.03903 X 10-18 

0 

0 

-3.4 

1.0 

-1.20061 X 10-19 

1.06197 X 10-19 

-1.14535 X 10-17 

0 

0 

-3.2 

1.0 

4.48104 X 10-19 

-3.96361 X 10-19 

-1.26253 X 10-18 

0 

0 

-3.0 

1.0 

-1,672460 X 10-9 

1.479340 X 10-9 

-1.39171 X 10-15 

0 

0 

-2.8 

1,0 

6.242140 X 10-9 

-5.521350 X 10-9 

-1.53411 X 10-14 

0 

0 

-2.6 

1.0 

-2.329760 X 10-8 

2.060740 X 10-8 

-1,69108 X 10-13 

0 

0 

-2.4 

1.0 

8.695370 X 10-8 

-7,691310 X 10-8 

-1.86410 X 10-12 

0 

-1.00000 X 10-11 

-2.2 

1.0 

-3.245380 X 10-'^ 

2.869740 X 10-'^ 

-2.05483 X 10-11 

-5.00000 X 10-11 

-4.25000 X 10-11 

-2.0 

1.0 

’1.211040 X 10-8 

-1.072200 X 10-6 

-2,26508 X 10-10 

-5.15000 X 10-10 

-5.15000 X 10-18 

-1.8 

1.0 

-4.523020 X 10-6 

3.989670 X 10-6 

-2,496840 X 10-9 

-5.685000 X 10-9 

-5.685000 X 10-9 

-1.6 

.9999999950 

1.684800 X 10-5 

-1.502460 X 10-5 

-2.752310 X 10-8 

-6.268750 X 10-8 

-6.268750 X 10-8 

-1.4 

,9999999490 

-6.324910 X 10-5 

5.460020 X 10-5 

-3.033920 X 10-7 

-6.910100 X 10-7 

-6.910100 X 10-7 

-1.2 

,9999994430 

2.320170 X 10-4 

-2.200570 X 10-4 

-3.344330 X 10-® 

-7.617000 X 10-8 

-7.616960 X 10“8 

-1.0 

.9999938560 

-9.105680 X 10-4 

6.419600 X 10-4 

-3.686440 X 10-5 

-8.395740 X 10-5 

-8.395270 X 10*8 

-.8 

.9999322760 

2.906350 X 10-3 

-4.371200 X 10-3 

-4.062620 X 10-4 

-9.247311) X 10-4 

-9.241600 X 10-4 

-.6 

.9992539710 

-.01624217000 

-5.226220 X 10-3 

-4.466160 X 10-3 

-.01010346100 

-.01003528000 

-.4 

.9918374290 

1.425470 X 10-3 

-.19277486800 

-.04778265700 

-.10135130100 

-.094 3 8690800 

-.2 

.9168273040 

-.50665763800 

-.33255690200 

-.38142198800 

-.51252817300 

-.30238007100 

0 

.5 

0 

.16384511900 

0 

0 

0 

.2 

.0831726960 

.50665763800 

.46984169900 

.38142198800 

.51252817300 

.30238007100 

.4 

8.16257 X 10-3 

-1.425470 X 10-3 

.02632517100 

.04778265700 

.10135130100 

.09438690800 

.6 

7,46029 X 10-4 

.01624217000 

9.038650 X 10-3 

4.466160 X 10-3 

.01010346100 

.01003528000 

.8 

6,77241 X 10-5 

-2.906350 X 10-3 

-9.405780 X 10-4 

4.062620 X 10-4 

9.247310 X 10-4 

9.241600 X 10-4 

1.0 

6.14417 X 10-6 

9.105680 X 10-4 

3.865550 X 10-4 

3.686440 X 10-5 

8.395740 X 10-5 

8.395270 X 10-8 

1.2 

5.57393 X 10-'^ 

-2.320170 X 10-4 

-9.136280 x 10-5 

3.344330 X 10-8 

7.617000 X 10-8 

7.616960 X 10-6 

1.4 

5.05660 X 10-8 

6.324910 X 10-5 

2.558630 x 10-5 

3.033920 X 10-7 

6.910100 X 10-7 

6.910100 X 10-7 

1.6 

4.58500 X 10-9 

-1.684800 X 10-5 

-6.754910 X 10-8 

2.752310 X 10-8 

6.268750 X 10-8 

6.268750 X 10-8 

1,8 

4.1600 X 10-19 

4.523020 X 10-8 

1.818960 x 10-8 

2.496840 X 10-9 

5.685000 X 10-9 

5.685000 X 10-9 ^ 

2,0 

3.7000 x 10-11 1 

-1,211050 X 10-8 

-4.865290 X IQ-^ 

2.26508 X 10-10 

5.15000 X 10-18 

5.15000 X 10-16 

2.2 

4.0000 X 10-12 

3,245490 X 10*’^ 

1.304300 X 10-7 

2.054 83 X 10-11 

4.62500 X 10-11 

4.62500 X 10-11 

2,4 

0 1 

-8,694880 X 10-8 

-3.493800 X 10-8 

1.86410 X 10-12 

5.00000 X 10-12 

5.00000 X 10-12 

2.6 

0 

2.329630 X 10-8 

9.360980 X 10-9 

1.69108 X 10-13 

0 

0 

2.8 

0 

-6.241790 X 10-9 

-2.508100 x 10-9 

1.53411 X 10-14 

0 

0 

3,0 

0 

1.672370 X 10-9 

6,71996 X 10-10 

1.39171 X 10-15 

0 

0 

3.2 

0 

-4.48079 X 10-19 

-1,80048 X 10-10 

1.26253 X 10-18 

0 

0 

3.4 

0 

1.20054 X 10-19 

4.82405 x 10-11 

1.14535 X 10-17 

0 

0 

3.6 

0 

-3,21662 X 10-11 

-1.29251 X 10-11 

1.03903 X 10-18 

0 

0 

3.8 

0 

8.61832 X 10-12 

3.46304 X 10-12 

9.42590 X 10-20 

0 

0 

4.0 

0 

-2.30912 X 10-12 

-9.27859 X 10-13 

8.55098 X 10-21 

0 

0 

4.2 

0 

6.18699 X 10-13 

2.48616 X 10-13 

7.75728 X 10-22 

0 

0 

4.4 

0 

-1,65825 X 10-13 

-6,66621 X 10-14 

7.03724 X 10-23 

0 

0 

4.6 

0 

4.46422 X 10-14 

1.80488 X 10-14 

6.38404 X 10-24 

0 

0 

4.8 

0 

-1.27539 X 10-14 

-5.76497 X 10-15 

5.79147 X 10-25 

0 

0 

5,0 

0 

0 

0 

5.25391 X 10 '26 
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36.000 5.273 x 10-3 5.239 x 10-3 6.307 x 10-3 6.264 x 10-3 3.316 x 10-2 3.293 x 10-2 3.601 x 10-1 3.571 x 10- 




TABLE 28.- COMPARISON OF RESULTS FOR THE SQUARE CAVITY 

FOR R = 10 


(a) Vorticity at center of moving wall 


Calculation method 

Points 

Vorticity at center 
of moving wall 

Spline 

15 X 15 

5.8884 

Finite difference 

15 X 15 

5.9264 

Finite difference, divergence form 

15 X 15 

5.9129 


(b) Maximum stream function 


Calculation method 

Points 

Maximum 
stream function 

Spline 

15 X 15 

-0.10027 

Finite difference 

15 X 15 

-.09790 

Finite difference, divergence form 

15 X 15 

-.09805 














TABLE 30.- COMPARISON OF RESULTS FOR THE SQUARE CAVITY FOR R = 100 






TABLE 31 .- COMPARISON OF RESULTS FOR THE VELCXITY u THROUGH POINT OF MAXIMUM STREAM FUNCTION FOR R = 100 
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TABLE 33.- COMPARISON OF RESULTS FOR THE 2x1 RECTANGULAR 

CAVITY FOR R = 100 

(a) Vorticity at center of moving wall 


Calculation method 

Points 

Vorticity 

^line 

Finite difference, divergence form 

29 X 15 
33 X 17 

7.1603 

7.3929 


(b) Upper vortex maximum stream function 


Calculation method 

Points 

Upper -vortex 
maximum 
stream function 

Spline 

29 X 15 

-0.10625 

Finite difference, divergence form 

33 X 17 

-.99286 

Reference 8 

21 X 21 

-.10204 


(c) Lower vortex maximum stream function 




Lower -vortex 

Calculation method 

Points 

maximum 
stream function 

^line 

29 X 15 

0.00094 

Finite difference, divergence form 

33 X 17 

.00059 

Reference 8 

21 X 21 

.00062 
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Figure 1,- Comparison of implicit spline and nondivergence finite -difference solutions 
with the exact solution to Burgers' equation for 51 points, = 1/24, equal spacing, 
and a = 0. 



Figure 2.- Comparison of two-step spline solution with the exact solution to Burgers' 
equation for 51 points, = 1/24, equal spacing, and a = 0, 
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' Exact solution 
O 2-Step spline solution 


Figure 3,- Comparison of two-step spline solution with the exact solution to Burgers 
equation for 51 points, u = 1/24, equal spacing, and a = 0. 



Exact solution 
Implicit spline solution 


Figure 4.- Comparison of implicit spline solution with the exact solution to Burgers 
equation for 51 points, u = 1/24, equal spacing, and a = 5. 




Figure 5.- Comparison of implicit spline solution with the exact solution to Burgers' 
equation for 37 points, p = 1/24, unequal spacing - more points near boundaries, 
and a = 0. 



Figure 6.- Comparison of implicit spline solution with the exact solution to Burgers' 
equation for 37 points, p = 1/24, unequal spacing - more points in corner region, 
and a = 0. 
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Exact solution 



Figure 7.- Comparison of implicit spline solution with the exact solution to Burgers 
equation for 31 points, u = 1/24, unequal spacing = 1.5, and ct = 0. 


Exact solution 


^ 2-Step spline solution 



-.6 -.4 


Figure 8.- Comparison of two-step spline solution with the exact solution to Burgers' 
equation for 31 points, v = 1/24, unequal spacing oj = 1.5, and ct = 0. 



Figure 9.- Comparison of two-step spline solution with the exact solution to Burgers' 
equation for 19 points, v = 1/24, unequal spacing Oj = 1.5, and a = 0. 



Figure 10.- Comparison of two-step spline solution with the exact solution to Burgers' 
equation for 19 points, u = 1/24, unequal spacing aj = 2, and o = 5. 
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Exact solution 


Implicit spline solution 


.0 -,8 -.6 -.4 -.2 


Figure 11.- Comparison of implicit spline solution with the exact solution to Burgers 
equation for 19 points, v = 1/24, unequal spacing = 1.75, and ct. = 5. 



— Exact solution 
O Implicit spline solution 


Small oscillations near boundary 


Figure 12.- Comparison of implicit spline solution with the exact solution to Burgers 
equation for 19 points, v = 1/24, unequal spacing = 2, and cr = 5. 





1 

Figure 13.- Comparison of implicit spline solution with the exact solution to Burgers' 
equation for 15 points, v = 1/24, unequal spacing = 1.5, and a = 4.5. 



Figure 14.- Comparison of implicit spline solution with the exact solution to Burgers' 
equation for 15 points, v = 1/24, unequal spacing m = 1.5, and cr = 7.5. 





Figure 15.- Comparison of two-step spline solution with the exact solution to Burgers 
equation for 15 points, v = 1/24, unequal spacing Oj = 1.75, cr = 0, and -3 ^ r/ ^ 3 



Figure 16.- Comparison of implicit spline solution with the exact solution to Burgers 
equation for 15 points, v = 1/24, unequal spacing aj = 1.75, and cr = 5. 
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Figure 17.- Comparison of implicit spline solution with the exact solution to Burgers’ 
equation for 15 points, w = 1/24, unequal spacing ffj = 2, and a = 5. 



Figure 18.- Comparison of the SADI solution with the exact solution to the two- 
dimensional diffusion equation. R = 1000; At = 9 x 10”^; 17 x 17 grid; 
unequal spacing; 0 ^ Y,Z ^ 4. 
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Nondivergence Form 

Finite Difference, 15 x 15 points 

Spline, 29 x 29 points 

Nondivergence Form 

Finite Difference, 57 x 57 points 

Divergence Form 

Finite Difference, 65 x 65 points 




t .4 .5 .6 -7 ,8 .9 1.0 


Figure 21.- Comparison of calculated velocity u through point of maximum 

for R = 100. 


Figure 22.- Comparison of calculated velocity u throv^h upper point of maximum 4 / 

for R = 100 and 2x1 rectangular cavity. 
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